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Abstract 

Background: The liver is central to most economically important metabolic processes in cattle. However, the 
changes in expression of genes that drive these processes remain incompletely characterised. RNA-seq is the new 
gold standard for whole transcriptome analysis but so far there are no reports of its application to analysis of 
differential gene expression in cattle liver. We used RNA-seq to study differences in expression profiles of hepatic 
genes and their associated pathways in individual cattle in either mild negative energy balance (MNEB) or severe 
negative energy balance (SNEB). NEB is an imbalance between energy intake and energy requirements for lactation 
and body maintenance. This aberrant metabolic state affects high-yielding dairy cows after calving and is of 
considerable economic importance because of its negative impact on fertility and health in dairy herds. Analysis of 
changes in hepatic gene expression in SNEB animals will increase our understanding of NEB and contribute to the 
development of strategies to circumvent it. 

Results: RNA-seq analysis was carried out on total RNA from liver from early post partum Holstein Friesian cows in 
MNEB (n = 5) and SNEB (n = 6). 12,833 genes were deemed to be expressed (>4 reads per gene per animal), 413 of 
which were shown to be statistically significantly differentially expressed (SDE) at a false discovery rate (FDR) of 
0.1% and 200 of which were SDE (FDR of 0.1%) with a >2-fold change between MNEB and SNEB animals. GOseq/ 
KEGG pathway analysis showed that SDE genes with >2- fold change were associated [P <0.05) with 9 KEGG 
pathways. Seven of these pathways were related to fatty acid metabolism and unexpectedly included 'Steroid 
hormone biosynthesis', a process which mainly occurs in the reproductive organs rather than the liver. 

Conclusions: RNA-seq analysis showed that the major changes at the level of transcription in the liver of SNEB 
cows were related to fat metabolism. 'Steroid hormone biosynthesis', a process that normally occurs in reproductive 
tissue, was significantly associated with changes in gene expression in the liver of SNEB cows. Changes in gene 
expression were found in this pathway that have not been previously been identified in SNEB cows. 
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Background 

Selective breeding for high milk yield in dairy cows has 
led to breeds in which the nutritional demands of the 
very high lactation rates following calving are in excess 
of that which the animal can metabolise from ingested 
feed [1]. This aberrant physiological state is known as 
negative energy balance (NEB) and is of particular 
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economic importance in the dairy industry because of 
its negative impact on health and fertility [1]. 

The principal physiological response of animals to 
NEB is to try and maintain homeostasis by mobilising 
body fat (and protein) reserves. This results in the re- 
lease of non-esterified fatty acids (NEFAs) from adipose. 
NEFAs are transported primarily to the liver where they 
are either fully oxidized to C0 2 , or converted to ketone 
bodies such as beta-hydroxybutyrate (BHB), or esterified 
into triacylglycerides (TAG) either for delivery into 
blood as very low density lipoproteins proteins (VLDL) 
or storage as cytosolic lipid droplets [2]. As with all 
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ruminants, dairy cattle have very low rates of VLDL syn- 
thesis and secretion [3] and are therefore susceptible to 
accumulation of high levels of TAG in the liver cells (re- 
ferred to as fatty liver or lipidosis) and high concentra- 
tions of ketone bodies in the blood (ketosis) [4] both of 
which are potentially detrimental to the health of the 
animals [5]. Although there are strategies to prevent and 
treat NEB, such as management of diet and inclusion of 
dietary supplements [6], [7] NEB and ketosis continue to 
have a significant negative economic impact in dairy cat- 
tle. This has led to an effort to understand NEB at the 
level of genes and their expression [8] . 

As with the majority of metabolic states, changes in 
gene expression in the liver are central to NEB. There 
have been several studies describing such changes at the 
level of transcription in the liver of NEB cows which 
have employed qPCR [9] and bovine microarrays 
[10,11]. However, microarray analyses rely on existing 
genome annotations and require design of a specialised 
array of probes based on information obtained from 
other methods such as sequencing in order to detect 
more complex regulations in gene expression such as al- 
ternative splicing [12]. Furthermore, the insensitivity of 
microarray platforms in detecting differential expression 
of certain metabolically important genes has been previ- 
ously highlighted by our own group who found that 
qPCR was able to detect changes in hepatic gene expres- 
sion SNEB cows which microarray analysis failed to de- 
tect [11]. This emphasises the need for a more sensitive 
technique for whole transcriptome analysis. 

RNA-seq is a relatively new technique that can be used 
to analyse changes in gene expression across the entire 
transcriptome [13,14], and is now being applied to a rap- 
idly increasing number of organisms [12]. This technol- 
ogy has distinct advantages over microarrays including 
the sensitive detection of all expressed genes without the 
need to generate an array of probes based on known se- 
quence, virtually no background noise and a much 
higher dynamic range. It has also led to the discovery of 
previously unidentified transcripts such as the recent 
discovery of enhancer RNA, a novel RNA species [15] 
and as the transcriptome is being sequenced, alternative 
splicing [16], single nucleotide polymorphisms (SNPs) 
[17], transcriptional fusions and chimeric transcripts 
[18] can be detected without the need to design an array 
of probes. 

Unlike microarrays, the bioinformatic tools for the 
analysis RNA-seq data are still at the early stages of de- 
velopment [19,20]. Partly for this reason, there are few 
reports to date on the use of RNA-seq for identification 
of significantly differentially expressed (SDE) genes be- 
tween physiologically different groups of animals and, so 
far, no such reports for cattle. The use of RNA-seq to 
identify differential gene expression in two pools of 



cattle embryos that were divergent for viability was re- 
cently reported [12]. However, as RNA from two pools 
of 20 embryos was analysed rather than individual 
embryos, the authors were unable to apply statistical 
analysis to their data to determine statistical significance 
of the differential gene expression. 

Here we describe the use of DEseq [21] and GOseq 
[22] to identify SDE genes and associated over-repre- 
sented biological KEGG pathways across the whole liver 
transcriptome of individual cattle. For this we compared 
RNA samples from liver of mild and severe NEB (MNEB 
and SNEB) cows. The objective of our study was to use 
RNA-seq to identify novel genes and their associated 
biological pathways which are important in SNEB in 
cattle. 

Results and discussion 

Library preparation and sequencing of polyA mRNA-seq 
libraries 

Twelve polyA RNA-seq single read libraries that had 
been prepared from total RNA extracted from liver from 
6 MNEB animals (libraries 1-6) and 6 SNEB animals (li- 
braries 7-12) were run as 36 bp single reads on 21 lanes 
randomly distributed across 3 flowcells (Table 1). Several 
libraries were run on one or more lanes at the same or 
different concentrations to determine the effect of li- 
brary concentration on cluster generation and to assess 
the consistency of sequence between different lanes and 
flow cells. For the 21 lanes, the average number of raw 
reads per lane was 13.66 million. Reads were aligned to 
the BCM4 genome assembly [23], using the ultrafast 
short read aligner, Bowtie (version 0.12.3) [24], Following 

Table 1 Summary of RNA-seq libraries and flow cells 

Lane Flow cell A Lane Flow cell B Lane Flow cell C 

1A library 1 1B library 12 1C library 2 

(4.5 pM) (5.5 pM) (5.5 pM) 

2A libraryl 2B library 6 2 C library 3 

(5.0 pM) (5.5 pM) (5.5 pM) 

3A library 5 3B library 2 3 C library 6 

(5.5 pM) (5.5 pM) (5.5 pM) 

4A Iibrary5 4B library 7 4 C library 7 

(6.0 pM) (5.5 pM) (5.5 pM) 

5A library 10 5B library 9 5C PhiX 

(4.0 pM) (5.5 pM) 

6A PhiX 6B PhiX 6C library 8 

(5.5 pM) 

7A library 10 7B library 11 7 C library 1 1 
(5.5 pM) (5.5 pM) (5.5 pM) 

8A library 10 8B library 4 8C library 12 



Libraries 1-6 = mild negative energy balance animals. Libraries 7-1 2 = severe 
negative energy balance animals. Libraries on the first flowcell were run at 
different concentrations (4.5-6.5 pM) to determine which concentration 
yielded optimum cluster yields. 
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removal of reads that either failed to align to BCM4 or 
aligned to more than one position in the genome, and 
the removal of all but one of the reads that mapped to 
exactly the same position on the genome (i.e. putative 
PCR duplicates), an average of 3.8 million reads per lane 
were retained for read abundance calculation. An over- 
view of these data is given in Additional file 1 and they 
have been deposited in NCBI's Gene Expression Omni- 
bus and are accessible through GEO Series accession 
number GSE37544 (http://www.ncbi.nlm.nih.gov/geo/ 
query/acc.cgi?acc=GSE37544). 

Transcriptional profile of cattle liver 

The lowest limit of detection was set to 5 or more 
uniquely aligned reads in at least one animal in either 
the SNEB or MNEB group. At this limit, 12,523 and 
12,833 genes were detected as expressed in MNEB and 
SNEB animals respectively. This is similar to the number 
of human Ensembl genes detected in human (12,191) 
and mouse (11,201) liver tissue by polyA RNA-seq [25]. 
RNA-seq analysis of mouse and human liver by other 
groups showed that the 10 most highly expressed genes 
make up 20-40% of the mRNA pool whereas other tissue 
transcriptomes e.g. brain, kidney and testis are more 
complex with the 10 most highly expressed genes con- 
tributing to just 5-10% of the total mRNA [25]. We 
found that in the liver of cows in NEB, on average 16% 
of the total RNA-seq reads aligned to only 10 genes and 
30% of the RNA-seq reads mapped to just 619 genes 
(Figure 1A) per lane. However, when we normalised for 
gene length by calculating values for fragments per kilo- 
base of exon per million fragments mapped (FPKM) to 
get a better indication of relative transcript number, the 
top 10 most highly expressed genes represented only 
0.77% of the transcript pool. There were highly abundant 
transcripts (FPKM >50) for approximately 1,200 genes 
in hepatic tissue analysed in the current study (Fig- 
ure IB). Depletion of very highly expressed long tran- 
scripts may be necessary to increase the sequencing 
coverage of less abundant shorter transcripts in mam- 
malian liver RNA-seq libraries. 

Lane effect 

To further determine whether the differential expression 
observed between MNEB and SNEB was due to a real 
difference between the groups and not simply random 
variation caused by running samples on different lanes 
(i.e. lane effect), correlation analysis was performed be- 
tween FPKM values from the two treatment groups and 
also the same 9 libraries run as technical replicates on 
two different lanes (Figure 2). This showed that the dif- 
ference in FPKM values (r 2 = 0.9783 NEB vs. SNEB 
r 2 = 0.9905 Lane A vs. Lane B) was predominately due to 
effect of treatment and not lane effect. 
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Figure 1 Complexity of NEB cattle liver transcriptome. (A) 

Average number of reads across 21 flowcell lanes aligning to genes. 

(B) Average FPKM values for genes across 21 flowcell lanes. 
I . / 

Identification of SDE genes 

DEseq is an R-based software package that was devel- 
oped specifically for the identification of SDE genes from 
raw counts of sequence reads generated by RNA-seq 
analysis that uniquely align to genes. RNA-seq library 1 
(corresponding to MNEB animal 1) was not included in 
these analyses as this animal subsequently became ill 
and a possible underlying infection may have skewed the 
gene expression data generated. DEseq analysis showed 
that 413 genes were SDE (FDR of 0.1%), 201 of which 
had a fold change of >2 (either up or down) between 
MNEB and SNEB animals (Additional file 2). 

Comparison of microarray and RNA-seq 

We previously analysed the same MNEB and SNEB liver 
RNA samples using the Affymetrix 23 K bovine gene 
microarray platform and, of the 5,229 genes which were 
detected as expressed, 416 were identified as SDE by the 
puma method using the recommended cut-off P-like 
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FPKMSNEB FPKMLaneA 

Figure 2 Lane effect. (A) Average FPKM values for each gene in SNEB animals correlated with average FPKM values for each gene in MNEB 
animals. (B) Average FPKM values for each gene for libraries 1, 2, 5, 6, 7, 10, 1 1, 12 run on lanes 1 A, 3B, 3A, 2B, 4B, 7A, 7B, 1B respectively 
correlated with average FPKM values for each gene for libraries 1, 2, 5, 6, 7, 10, 1 1, 12 run on lanes 2A, 1C, 4A, 3C, 4C, 8A, 5A, 7B, 8C respectively. 



value of P <0.05 [11]. Of these, only 56 had fold changes 
of >2 (up or down) in SNEB animals. 55 of the 416 
microarray and 413 RNA-seq SDE genes, and 27 of the 
201 RNA-seq and 56 microarray >2-fold SDE genes, 
were detected as SDE by both platforms (Additional file 
3). All genes detected as SDE on both platforms had fold 
changes in the same direction. This approximate com- 
parison of the two data sets shows that, particularly for 
the lower fold changes, many genes were detected as 
SDE on one platform but not the other. However, RNA- 
seq detected expression of more than twice as many 
genes as the 23 K bovine microarray, and nearly four 
times as many >2-fold SDE genes. 

Physiology of SNEB model 

Soon after calving in high-yielding dairy cows, nutri- 
tional and energetic demands can increase 4-fold within 
a single day as the animal undergoes a physiological up- 
heaval to meet the demand by the mammary gland for 
substrates, particularly glucose, for milk production. The 
main changes of this upheaval are a reduction of insulin 
concentration in the blood, a reduction of lipogenesis in 
adipocytes, an increase in export of fats from adipose to 
liver, and a dramatic increase in gluconeogenesis in liver. 
Unlike non-ruminants, where glucose is absorbed via 
the small intestine, approximately 90% of all glucose in 
ruminants is synthesised by gluconeogenesis in the liver 
from volatile fatty acids (VFAs) derived from bacterial 
fermentation in the rumen. VFAs are absorbed across 
the rumen wall and into the hepatic portal vein which 
delivers them to the liver. The increase in hepatic gluco- 
neogeneis at the onset of lactation leads to production 
of very high quantities of glucose. The glucose entry into 
the blood stream in top performing cows producing 
90 kg milk per day has been estimated at 7.4 kg per day 
of which 4.4 kg will end up as lactose in milk [26] . 

In our model we were comparing two groups of lactat- 
ing animals in the early post-calving period (13 or 



14 days post-calving). One of these groups was in MNEB 
and the other was in a state of SNEB which was artifi- 
cially induced by restricting feed and trebling the milk- 
ing frequency. From this we hoped to identify novel 
genes and pathways which are important in SNEB in 
high-yielding dairy cows. Hormones, metabolites and 
aspects of liver composition measured at the time of tis- 
sue collection showed that, as predicted, relative to 
MNEB cows, SNEB cows had undergone significant 
(P <0.05) changes in energy balance (-3.6-fold), blood 
glucose (-1.5-fold), blood NEFA (+ 4.7-fold), blood BHB 
(+7.4), blood IGF1 (-4.8-fold), liver TAG (+3.9-fold), liver 
glycogen (-10.4-fold). Blood insulin, oestradiol, urea and 
liver weight were not different between SNEB and 
MNEB animals [9]. The level of cholesterol was also 
lower in the SNEB cows (Additional file 4). These 
changes are widely known to occur in SNEB cows [10] 
where there is an increase in the export of NEFAs from 
adipose to the liver which exceeds the normal oxidative 
capacity of the liver. This results in an increase in ketone 
bodies (BHB) in the blood resulting from partially oxi- 
dised NEFAs, and an accumulation of TAGs in the liver. 
In addition, IGF1 levels drop as hepatocytes become re- 
fractory to growth hormone (GH). Under normal condi- 
tions, growth hormone (GH) is released from the 
anterior pituitary in the brain, binds to the growth hor- 
mone receptor (GHR) on the surface of hepatocytes, and 
this binding triggers a signalling cascade that leads to 
increased IGF1 synthesis by the liver for export to the 
blood. In NEB cows the liver becomes refractory to GH 
and this results in low circulating IGF1 [10]. 

Over-represented KEGG pathways 

In our RNA-seq study, nine KEGG pathways were asso- 
ciated with >2-fold SDE genes most of these were 
related to metabolism of fats (Table 2, Figures 3, 4, Add- 
itional file 5). Pathways related to the metabolism of car- 
bohydrates were only over-represented when SDE genes 
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Table 2 Summary of KEGG pathways associated with >2 
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Fat digestion and absorption 


0.004934 


02010 


ABC transporters 
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a-linolenic acid metabolism 
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Biosynthesis of unsaturated 
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Pathway maps for KEGG pathways associated with >2-fold SDE genes 
(FDR0.1%) that are not included in the main text are in Additional file 5. 



with < 2-fold change were included in the pathway ana- 
lysis, (Table 3, Additional file 6). This suggests that in 
our model, at the level of transcription in the liver, to 
compensate for the reduced supply of VFAs from the 
liver due to feed restriction, the major changes in hep- 
atic gene expression in SNEB cows were related to the 
alternative supply of carbon from fat for gluconeogen- 
esis. The most interesting finding from our RNA-seq 
analysis was that 'Steroid hormone biosynthesis', a 
process which is normally exclusive to the gonads and 
adrenal glands, was the KEGG pathway most signifi- 
cantly associated with >2-fold SDE genes in SNEB cow 
liver. Details of this and other selected pathways are dis- 
cussed below in the context of the expected physio- 
logical changes in the liver of SNEB cows. 

Carbohydrate metabolism 

The 1.5-fold drop in blood glucose in SNEB cows would 
suggest decreased gluconeogenesis in the liver of these 
animals. However, only a minor decrease (<2-fold) of 
pyruvate carboxylase (PYK (KEGG pathway map) = 
PKLR (SDE gene list), one of the major rate limiting 
enzymes in gluconeogenesis [27] was identified in the 
KEGG Insulin signalling pathway (Additional files 2 and 
6). The KEGG pathway 'Glycolysis/gluconeogeneis' was 
not identified as over-represented. This indicates that, in 
feed-restricted SNEB animals, there was only a slight de- 
crease in the rate of gluoneogenesis and that demands 
for carbon for this process were mostly being met by 
oxidation of NEFAs transported from adipose. There 
were no gene pathways to related to the 10-fold drop in 
liver glycogen in observed in SNEB cows even though 
genes coding for enzymes involved in glycogen metabol- 
ism (e.g. PYG, PHK) were detected as expressed in the 
KEGG Insulin signalling pathway (Additional file 6). 



Major changes in transcription of genes controlling gly- 
cogenosis may have occurred prior to collection of liver 
tissue for RNA extraction. 

Uncoupling of the GH/IGF1 axis 

Although pathways related to the GH/IGF1 axis were 
not identified as overrepresented among the SDE genes 
identified by RNA-seq, several genes within in this path- 
way were detected as SDE (GHR -1.6-fold, IGFALS -5- 
fold, IGF1 -5-fold) (Additional file 2). The KEGG p53 
signalling pathway', which is related to stress and 
includes IGF1, was overrepresented (Table 3). Hepatic 
expression of IGF1, GHR and IGFALS is well known to 
be reduced in SNEB cows and has been previously 
described [9,11]. Interestingly, for the same total RNA 
samples, IGF1 was detected as SDE by both RNA-seq 
and qPCR but not by microarray (Additional file 3) [11], 
even though the amount of blood IGF1, which is made 
mainly in the liver, was significantly reduced in SNEB 
cows. Microarray also failed to detect GHR as SDE. 
However, our microarray study did, whereas our RNA- 
seq study did not, identify IGFBP2 and IGFBP3 (genes 
encoding IGF1 binding proteins) as SDE (<2-fold) in 
SNEB cow liver. 

TAG accumulation and increased blood NEFA 

Increased TAG accumulation in the liver was reflected by 
overrepresentation of the KEGG 'PPAR signalling path- 
way (Table 2, Figure 3) as angiopoietin-like 4 gene (PGAR 
(KEGG pathway map, Figure 3)=ANGPTL4 (SDE gene 
list, Additional file 2) was upregulated in this pathway. 
Angiopoietin-like protein 4 has emerged as a key regulator 
of plasma cholesterol, triglyceride, and NEFA concentra- 
tions and, in cattle, is most highly expressed in the liver 
and adipose [28]. The upregulation of ANGPTL4 due to 
fasting, corresponding to the SNEB condition, was also 
reported in adipose tissues of lactating goats [29]. Upregu- 
lation of ANGPTL4, which is stimulated by PPAR ago- 
nists, leads to inhibition of lipoprotein lipase activity in 
adipose, reduced VLDL-TAG utilisation, and increased 
lipolysis. This gene was previously identified by microarray 
to be up-regulated in the liver of induced ketosis cows 
[10] and SNEB cows [11] and it was proposed that 
increased ANGPTL4 serves as a signal for lipolysis and 
contributes to sustained release of NEFA and lipid accu- 
mulation in the liver [10]. 

Oxidation of NEFAs and ketogenesis 

Due to the higher amount of NEFAs in the blood of 
SNEB animals, we were expecting changes in expression 
of genes to reflect increased oxidation of fatty acids in 
SNEB animals. However, although eight of the nine fatty 
acid oxidation genes in the KEGG PPAR signalling path- 
way were detected, none were SDE (Additional file 2). 
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Figure 3 Changes in hepatic gene expression in the KEGG PPAR signalling pathway in SNEB compared to MNEB cows. Blue boxes are 
genes detected by RNA-seq as expressed but not SDE, red boxes are SDE genes with reduced expression in SNEB animals, green boxes are SDE 
genes with increased expression in SNEB animals. 



Acetyl- Co A carboxylase beta (ACACB) which is thought 
to control fatty acid oxidation was SDE (upregulated 3- 
fold) in RNA-seq but not picked up in the pathway 
analysis. KEGG PPAR signalling pathway is also 
involved in ketogenesis but, although the main gene 
involved in this (HMGCS2) was detected, it was not 
identified as SDE. This may be because at the time of 
slaughter there was little difference in the rate of oxida- 
tion and ketogenesis between the two groups, or that 
there were differences but not at the level of transcrip- 
tion. OLR1 (oxidized low density lipoprotein (lectin-like) 
receptor 1) is included in the KEGG PPAR signalling 
pathway and was detected as upregulated (although 
read counts were very low for this gene). The protein 
encoded by OLR1 binds, internalizes and degrades oxi- 
dized low-density lipoprotein (oxLDL) and is thought to 



have a role in lowering systemic oxidative stress by re- 
moval of oxLDL from blood. 

Impaired PUFA synthesis 

The most SDE gene we identified by RNA-seq was FADS2 
(previously named A6 desaturase) which was down-regu- 
lated in SNEB animals (Figure 3, Additional files 2) and fea- 
tured in 3 of the nine >2-fold KEGG pathways (a-Linolenic 
acid metabolism, Biosynthesis of unsaturated fatty acids, 
and PPAR signalling) (Additional file 5). FADS1 was also 
shown to be down-regulated in SNEB animals (Additional 
files 2). FADS1 and FADS2 were previously identified by 
microarray as downregulated in liver of cows with induced 
clinical ketosis (days 9-14 post partum) [10]. However, in 
our microarray study, FADS2, but not FADS1, was identified 
as SDE. FADS2 catalyses the initial desaturation step of a- 
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Figure 4 Changes in hepatic gene expression in the KEGG Steroid hormone biosynthesis pathway in SNEB compared to MNEB cows. 

Blue boxes are KEGG enzymes coded for by genes detected by RNA-seq as expressed but not SDE, red boxes are KEGG enzymes coded for by 
SDE genes with reduced expression in SNEB animals, green boxes are KEGG enzymes coded for by SDE genes with increased expression in SNEB 
animals. KEGG enzyme 1.14.15.6 is coded for by HSA CYP11A1, KEGG enzyme 1.14.13.17 is coded for by HSA CYP7A1, KEGG enzyme 2.4.1.17 is 
coded for by HSA UG2A1, KEGG enzyme 2.8.2.4 is coded for by HSA SULT1E1. 
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Table 3 Summary of KEGG pathways associated with all 



SDE genes 


Pathway 
code 


Pathway name 


Over represented 
P value 


04976 


Bile secretion 


0.001483 


00140 


Steroid hormone biosynthesis 


0.002046 


04950 


Maturity onset diabetes 
of the young 


0.004947 


00910 


Nitrogen metabolism 


0.007765 


00061 


Fatty acid biosynthesis 


0.011248 


loz I y 


PPAR signalling pathway 


U.U I o I I z 


03320 


Bladder cancer 


0.018941 


04910 


p53 signalling pathway 


0.020756 


04115 


Insulin signalling pathway 


0.021653 


04977 


Vitamin digestion and absorption 


0.026059 


00592 


a-linolenic acid metabolism 


0.030405 


04975 


Fat digestion and absorption 


0.034363 


04080 


Neuroactive ligand-receptor 
interaction 


0.040178 



Pathway maps for KEGG pathways associated with all SDE genes (FDR 0.1.%) 
are in Additional file 6. 



linolenic acid (C18:3) and linoleic acid (C18:2) in the en- 
zyme cascade that ultimately leads, via FADS1, to synthesis 
of the long chain polyunsaturated fatty acids (LC-PUFAs) 
arachidonic acid (C20:4) and eicosapentaenoic acid (C20:5). 
Synthesis of LC-PUFAs arachidonic acid and eicosapenta- 
enoic acid mainly occurs in the liver and adipose [30] and 
suppression of their synthesis in SNEB cows may contribute 
to the poor fertility in high-yielding dairy herds. It was re- 
cently reported that polymorphisms in FADS1 and FADS2 
were associated with reduced serum arachidonic levels in 
humans [31]. Also, deletion of FADS2 gene expression in 
mice leads to sterility in both male and female mice without 
affecting their viability [32]. 

Unexpected overrepresentation of KEGG pathway 'steroid 
hormone biosynthesis' 

One of the most interesting findings of our RNA-seq 
study was that the KEGG pathway 'Steroid hormone bio- 
synthesis' (Table 2, Figure 4) showed the most significant 
association with > 2-fold SDE genes in SNEB animals. 
This is unexpected as, under normal conditions, steroid 
hormone biosynthesis occurs almost exclusively in the 
adrenal glands and gonads, whereas liver is the site for 
steroid hormone inactivation. Four genes were changed 
in this pathway. CYP11A1 was up-regulated while 
UGT2al, SULTE1 and CYP7A1 were down- regulated. 
CYP11A1 was also identified as SDE in SNEB cow liver 
in our microarray study [11] but not in the liver of 
induced ketosis cows [10], CYP7A1, UGT2A1 and 
SULT1E1 were not identified as SDE in either of those 



studies. CYP7A1 catalyses the rate limiting step of con- 
version of cholesterol to bile acids and was also included 
in the KEGG pathway 'Bile secretion' which was overre- 
presented in the > 2-fold SDE genes identified by RNA- 
seq (Table 2, Additional file 5). UGT2A1 and SULT1E1 
inactivate oestrogens by glucuronidation and sulphation 
respectively, prior to their excretion. The down-regula- 
tion of these genes in the liver of SNEB animals indi- 
cates that bile synthesis and inactivation of oestrogens in 
liver is reduced in SNEB possibly to conserve cholesterol 
which was reduced in SNEB animals (Additional file 4). 
This also suggests that lower blood concentration of oes- 
trodiol in preovulatory high yielding dairy cows may not 
be due to increased oestrogen inactivation in liver. 

CYP11A1, a novel gene in metabolic adaptations to post- 
partum bovine SNEB 

Although we detected upregulation of CYP11A1 in liver 
of SNEB cows in our microarray study, 'Steroid hormone 
biosynthesis' was not associated with the microarray 
SDE genes, probably because CYP7A1, UGT2A1 and 
SULTE1 were not identified as SDE, so the implications 
of upregulation of CYP11A1 in SNEB cows were not dis- 
cussed [11]. To the knowledge of the authors, other than 
in our own microarray study, expression of CYP11A1 
has not been previously described in adult cow liver. 
Other workers found CYP11A1 to be undetectable by 
qPCR in the liver of adult cows and consequently used 
liver as a negative control for expression studies of 
CYP11A1 in granulosa [33]. However, CYP11A1 does ap- 
pear to be expressed in foetal and juvenile cows [34]. 
Up-regulation of this gene is unexpected as it catalyses 
the initial step of conversion of cholesterol into pregne- 
nolone which is then converted by CYP17A1 to dehy- 
droepiandrosterone (DHEA). DHEA is the precursor for 
a number of steroid hormones including oestrogen. Al- 
though other tissues, such as brain and foetal rat liver, 
also express the enzymes for steroid biosynthesis [35] 
these enzymes are not normally expressed in adult mam- 
malian liver. Unexpected up-regulation of CYP11A1 was 
also observed in vitro in human hepatocytes after treat- 
ment with adenovirus expressing PGC-1 alpha [35] 
which is a key integrator of many of the signalling path- 
ways that are induced in the liver and muscle upon fast- 
ing. In that study, CYP17A1 was also upregulated and 
the authors proposed that DHEA (or one of its metabo- 
lites) in liver could be involved in a novel hepatic signal- 
ling pathway or may serve to protect hepatocytes from 
the increased oxidative state brought on by fasting. 
DHEA has been shown to protect endothelial cells from 
apoptosis. They also suggested that CYP11A1 expression 
in feed restricted hepatocytes may be involved in regula- 
tion of cholesterol homeostasis as CYP11A1 catalyses 
modifications on cholesterol to produce oxysterols 
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which are known ligands for liver X receptors (LXRs). 
LXRs regulate cholesterol homeostasis and lipid metab- 
olism. It would therefore be interesting to see if DHEA 
levels in SNEB cow liver are elevated. In both our micro- 
array and RNA-seq study, SDE genes included in the 
IPA canonical signalling pathway 'LXR/RXR activation' 
were changed (ABCG8, ACACA (downregulated); ARG, 
CCL2, IL1RN, RXRG (upregulated)). In our RNA-seq 
data, many of the genes, including CYP17A1, in the 
KEGG pathway 'Steroid hormone biosynthesis' were not 
detected as expressed (Additional file 5). This may either 
be due to the fact that these genes are duplicated in the 
genome and thus reads were eliminated because they 
mapped to more than one locus (e.g. CYP17A1 (human 
gene name) has an exact copy (LOC784256) at another 
position in the Bos taurus genome), transcript abun- 
dance was too low for detection, or that they are not 
expressed in the liver thus preventing complete synthesis 
of steroid hormones in this organ. 

Conclusions 

RNA-seq analysis showed that the major changes at the 
level of transcription in the liver of SNEB cows were 
related to fat metabolism. Unexpectedly, one of the most 
significantly changed biological pathways was 'Steroid 
hormone biosynthesis' where CYP11A1 was upregulated. 
In adult cows, this gene is normally only expressed in 
the adrenal glands and ovaries and not liver. CYP11A1 is 
a potentially novel gene that may play a key role in 
metabolic adaptations to NEB in high-yielding dairy 
cows. 

Materials and methods 

Animal model 

The animal model employed in this study has been 
described previously [9], and all procedures were carried 
out under license in accordance with the European 
Community Directive 86-609-EC. The nutritional and 
lactational management regime employed were designed 
to create significant divergence in the energy balance 
(EB) profiles of the cows in early lactation. In brief, mul- 
tiparous Holstein-Friesian cows (n = 24) were blocked 2 
wk prior to expected calving date according to parity, 
body condition score (BCS), and previous lactation yield 
(average lactation 6,477 ±354 kg) and randomly allo- 
cated to mild (MNEB, n = 12) or severe (SNEB, n = 12) 
NEB groups. MNEB cows were fed ad libitum grass sil- 
age and 8 kg/day concentrates and milked once daily; 
SNEB cows were fed 25 kg/day silage and 4 kg/day con- 
centrate and milked three times daily. Measurements of 
BCS and EB were used to select cows that showed 
extremes in EB from each group (MNEB, n = 5; SNEB, 
n = 6). Cows were slaughtered on days 6-7 of the first 
follicular wave after calving (mean number of days 



postpartum: MNEB mean 13.6 ±0.75, range 11-15; 
SNEB mean 14.3 ±0.56, range 13-16), based on daily 
transrectal ultrasonography. 

Collection of liver tissue for RNA and TAG analysis 

The entire liver was removed within 15-30 min after 
slaughter and weighed. Samples weighing approximately 
1 g were dissected, rinsed in RNase-free phosphate buf- 
fer, snap-frozen in liquid nitrogen, and stored at -80°C. 
For triacylglyceride (TAG) analysis, total lipids were 
extracted from 50 mg samples of liver as previously 
described [9]. 

Blood sampling and metabolite assays 

Stabilized (EDTA-treated) whole blood samples were 
collected on the day of slaughter by jugular venipuncture 
for haematological analysis. Blood samples were ana- 
lyzed for glucose, NEFA, (3-hydroxybutyrates (BHB), and 
urea using appropriate kits and an ABX Mira autoanaly- 
ser (ABX Mira, Cedex, France). All metabolite assay 
coefficients of variation were low and typically <5%. 

RNA extraction and quality analysis 

Total RNA was prepared from 100-200 mg of fragmen- 
ted frozen liver tissue using the TRIzol reagent (Sigma- 
Aldrich, Dorset, UK). Tissue samples were homogenized 
in 3 ml of TRIzol reagent and chloroform and subse- 
quently precipitated using isopropanol (Sigma Chemical, 
Wicklow, Ireland). RNA samples were stored at -80°C. 
20 ug of total RNA from each sample for genomic DNA 
contamination was treated with the RNase-free DNase 
set (QIAGEN, Crawley, West Sussex, UK) and purified it 
using the RNeasy mini kit in accordance with guidelines 
supplied (QIAGEN, Crawley, West Sussex, UK). RNA 
quality and quantity were assessed using automated ca- 
pillary gel electrophoresis on a Bioanalyzer 2100 with 
RNA 6000 Nano Labchips according to manufacturers 
instructions (Agilent Technologies Ireland, Dublin, Ire- 
land). Samples of RNA had 28 S/18 S ratios ranging 
from 1.8 to 2.0 and RNA integrity number values of be- 
tween 8.0 and 10.0. 

Preparation and sequencing of polyA mRNA-seq ilumina 
libraries 

polyA RNA was isolated from 5-10 ug total RNA with 
oligo(dT) beads using two rounds of oligo-dT purifica- 
tion. 5-10 ug RNA was fragmented with zinc fragmen- 
tase (Applied Biosystems, Warrington, UK), first strand 
cDNA synthesis was performed using the Invitrogen 
random hexamer primers and Superscript II (Invitro- 
gen), second strand synthesis was performed using Invi- 
trogen DNA Polymerase 1 (Invitrogen). End repair and 
polyadenylation were performed using NEB Next Tailing 
Module (New England Biolabs). ® End Repair Module 
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and NEB Next dA- Illumina single read adapters were 
ligated to blunt ended, polyadenylated fragments with a 
NEB Quick ligation kit (New England Biolabs). Adapter- 
ligated cDNA fragment libraries were run on an Illumina 
GAII using version 3 sequencing kits and version 3 sin- 
gle read cluster generation kits. 

Read alignment 

The reads from each of the lanes were aligned separately 
to the BCM4 genome assembly [23], using the ultrafast 
short read aligner, Bowtie (version 0.12.3) [24]. Fastq 
output files from the sequencer were used as input, spe- 
cifying the following options: quality scores are ASCII 
characters equal to the Phred quality scores plus 64 (- 
-Solexal.3-quals); the maximum number of mismatches 
allowed in the first 28 bases is 2 (-n 2, -1 28); suppres- 
sing all alignments for any read that has more then 1 re- 
portable alignment (-m 1); retained alignments were 
reported in SAM format (-S). 

Read abundance calculation 

The abundance of mRNAs for all annotated genes from 
the ENSEMBL v59 annotation of the bovine genome 
[36], was calculated using the software package HTseq 
(version 0.4.4p6) (http://sourceforge.net/projects/htseq/). 
The script HTseq-count provided with this package was 
used to count the number of reads that mapped to each 
annotated gene, allowing, in some cases, for reads to 
partially overlap with the exons and still be counted for 
that gene (i.e. -m union). As our sequencing analysis 
included technical replicates of individual samples, we 
summed the counts for all these lanes, resulting in a sin- 
gle count for each gene for each sample. The counts for 
all the samples were collated into one file and any gene 
with fewer than 5 reads in all samples was excluded 
from the subsequent statistical analysis of differential 
gene expression. 

Identification of SDE genes and pathway analysis 

We identified SDE using DEseq (version 1.1.11) [21]. 
DEseq uses a generalisation of the Poisson model, the 
negative binomial distribution, to model biological and 
technical variance and test for differential expression be- 
tween two experimental conditions. As there are nearly 
30,000 genes annotated in the bovine genome, the statis- 
tical tests in each analysis were corrected for multiple 
testing using the Benjamini and Hochberg (BH) method 
[37] as implemented in R (version 2.12.0). All genes that 
were found to be SDE between the two experimental 
conditions (at an FDR p-value cut-off of 0.1) were 
retained for further analysis. 

GOseq and KEGG pathway analysis tools were used to 
identify biological pathways that were significantly 
enriched in the data set of SDE genes (Table 2, 



Additional file 3). To facilitate this, reads were firstly 
converted to their human orthologs. GOseq specifically 
adjusts for the higher abundance of reads from long or 
highly expressed genes in RNA-seq experiments when 
assessing a gene list for over-representation of biological 
functions [22]. Genes were mapped to the KEGG data- 
base [38] for pathway analysis using GOseq [22]. 

Additional file 



Additional file 1. Table of number and assignment of reads for 
each flowcell lane. 

Additional file 2. Table SDE genes and their fold-change. 

Additional file 3. Comparison of microarray and RNA-seq SDE 
genes. 

Additional file 4. Graph of effect of SNEB compared to MNEB on 
blood concentrations of major NEB associated metabolites 1- 
14 days post-calving. 

Additional file 5. Maps of KEGG pathways (associated with >2 fold 
(FDR 0.1%) SDE genes) overrepresented in SNEB animals. 

Additional file 6. Maps of KEGG pathways associated with all SDE 
(FDR 0.1%) genes. 
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